Multiple liquid-liquid critical points and density anomaly in core-softened potentials 
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The relation between liquid-liquid phase transitions and waterlike density anomalies in core- 
softened potentials of fluids was investigated in an exactly solvable one dimensional lattice model 
and a in a three dimensional fluid with fermi-like potential, the latter by molecular dynamics. 
Both systems were shown to present three liquid phases, two liquid-liquid phase transitions closely 
connected to two distinct regions of anomalous density increase. We propose that an oscillatory 
behavior observed on the thermal expansion coefficient as a function of pressure can be used as a 
signature of the connection between liquid-liquid phase and density. 
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The phase behavior of single component systems as 
particles interacting via the o-called core-softened (CS) 
potentials has received attention since the pioneering 
work of Stell and Hemmer proposing the possibility of 
a second critical point in addition to the usual liquid- 
gas critical point*. These potentials exhibit a repul- 
sive core with a softening region with a shoulder or 
a rampir— which are analytically and computationally 
tractable while being capable of retaining the qualita- 
tive features of the real fluid systems. Furthermore, 
Debenedetti and collaborators 4 using thermodynamic ar- 
guments showed that the isobaric thermal expansion co- 
efficient (a) of these potentials might have an anomalous 
negative value and consequently a region where density 
increases with temperature. Since for high temperatures 
density decreases with temperature these potentials can 
exhibit a temperature of maximum density (TMD) con- 
necting the two regions. 

The thermodynamic anomalies predicted by these 
models occur in liquids such as Te*, Ga, Bi*, S*»*, liquid 
water—, and GeisTesg^, and were found in simulations 
for silic a 11 ' 12 , silicon^*, and BeF2**. In addition experi- 
ments in phosphorous indicate the presence of a liquid- 
liquid phase transitions and similar transitions were ob- 
served by simulations in models of silica**, silicon** and 
liquid water—. 

In the particular case of water the hypothesis of the ex- 
istence of two liquid phases has been indirectly supported 
by experimental results in confined systems**. In spite 
of the limit of 235 K below which water cannot be found 
in the liquid phase without crystallization, two amor- 
phous phases, a low density amorphous phase and a high 
density amorphous phase, were observed at much lower 
temperatures and their relation to a metastable liquid- 
liquid phase transition was argued**. More recently a 
third amorphous phase, the very high density amorphous 
phase, has been observed** what suggests the possibility 



of the existence of the very high density liquid phase. 

Therefore the issue of a liquid-liquid phase transition 
and its connection with the existence of a region in the 
pressure-temperature phase diagram where the density 
decreases with the decrease of temperature is itself an 
interesting topic. It is accepted that the presence of two 
accessible length scales in the potential allow for the sys- 
tem to have two liquid phases and a density anomaly. 
The accessibility is the ingredient that explains why a 
density anomaly derived in ID does not necessarily hold 
at higher dimensions and why a density anomaly derived 
for a smooth potential might be lost if the slope linking 
the two length scales becomes infinite** - — . 

In this letter we show, by means of an exactly solv- 
able ID model and numerical simulations of a similar 
3D potential, that a three length scales potential might 
exhibit three critical points and two density anomalous 
regions if the different length scales would be accessible. 
To the best of our knowledge this is the first time that 
core softened potentials are shown to have three critical 
points and two temperature of maximum density lines. 
In addition we propose propose a new way to identify if 
a liquid-liquid critical point would be present by explor- 
ing the behavior of a for temperatures above the critical 
temperature. We show that, for systems in which there is 
a density anomalous region close to criticality, a exhibits 
a peculiar behavior diverging to +oo for pressures above 
the critical pressure and to — oo for pressures below the 
critical pressure. 

The core-softened potentials we analyze are illustrated 
in Fig. Q] in units of length a and of energy e. They are 
composed by a hard core of diameter a, two repulsive 
shoulders and an attractive well with depth — Ae where 
we employ A = 1/10 for the ID case and A = 1 for the 
3D case. Both potential are detailed below. 

We first investigate the ID to obtain analytical insights 
and to become familiar with the properties that will also 
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FIG. 1. Representation of the pair potentials investigated on 
this paper: a one dimensional lattice fluid and an isotropic 
three dimensional continuous core softened model with three 
scales of interaction. Top (bottom) horizontal axis used to 
represent continuous (lattice) system. 



appear in the 3D case. The linear lattice has length L 
and N sites separated by distance a — L/N. The exact 
form of the Gibbs free energy derived in the framework 
of the Takahashi method^ is given by 



G(T,P,N) = -Nk B T In 



-Ph{r;P) 



x=l 



(1) 



where j3 = l/(fc^T) with ks the Boltzmann constant, 
and h(r; P) = U(r) + Pr the microscopic pair enthalpy, 
with U (r) being the interaction energy between neighbor 
molecules. Expression ([T]) is calculated using U (a) = 3e, 
U(2a) = e, U{3a) = -e/10, and U{r) = for r > 3cr. 

The minimization of the Eq. (JT|) in the ground state 
results in three configurations depicted in Fig. [TJ a 
low density liquid (LDL) phase, a high density liquid 
(HDL) phase and a very high density liquid (VHDL) 
phase, besides a gas phase (G), not illustrated. Co- 
existence between these phases allows for three ground 
state phase transitions (GSTP) to occur at reduced pres- 



sures (P* = Poje) P* G _ LDL 
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and P, 



HDL — V HDL 



2.0. In what follows temperature is 
reduced as T* = k B T/e. 

The full phase diagram for the ID case is illustrated on 
Fig. [5] (a) with symbols indicating the three GSPT de- 
scribed in the last paragraph. From the LDL-HDL crit- 
ical point (circle) emerges a line of temperature of mini- 
mum density (TmD) followed by a line of temperature of 
maximum density (TMD), while from the HDL- VHDL 
critical (triangle) point appears a TMD. These lines sep- 
arate the region where the thermal expansion coefficient 
is positive from the regions where a is negative. Fig. [5] 
(b) illustrates the density as function of the temperature 
for various pressures, some of them crossing the TMD 
and the TmD lines. 

Further comprehension on the relation between density 
anomaly and GSPT can be gained by visualizing a as a 




FIG. 2. (a) Pressure vs. temperature phase diagram for 
the ID case. Continuous and dotted lines indicate the TmD 
and the TMD lines, while symbols locate the ground state 
phase transitions, (b) Density as a function of temperature, 
for various pressures. The inset shows the reduced volume as 
a function of the reduced temperature near the two critical 
points, according to ©. 
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FIG. 3. Thermal expansion coefficient (a* — ea) as a function 
of pressure at fixed temperatures, for the ID case. 



function of pressure (at fixed temperature), as shown in 
Fig. [3J In usual systems a — > oo as the critical point is 
approached from any path in the pressure-temperature 
plane. The peculiarity of models in which the criticality 
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is associated with the density anomaly is that below the 
critical pressure P c the sign of a is negative, while for 
P > P c its sign is positive. Consequently a displays 
an oscillatory behavior as a function of pressure, with 
a(P — > P~) -» -oo and a(P — > P+) — > +00, which is 
a signature of the connection between the TMD line (or 
TmD) and the critical point^l. 

The relation between criticality and anomalous density 
behavior can be rationalized by looking at the thermo- 
dynamic close to the GSPT within the approximation 
adopted in Refi2£. For a lattice model with core softened 
interaction the reduced Gibbs free energy per particle in 
the vicinity of any GSPT can be written in terms r\ and 
T2 , the average distances between the particles in the two 
coexisting phases, leading to 



]{t,p) = h c +p - tin 



2 cosh ( — 
t 



(2) 



where g = G/(NP c r c ). 
and h c = h(r\)/(P c r t 
discussed above: r\ = 
n ) = 5, and r c = (n -+ 



t = k B T/P c r c , p={P- P c )/P c , 
). For the LDL-HDL transition 
2cr, r 2 = 3cr, d = (n + r 2 )/(r 2 - 

r 2 )/2 = 5cr/2. It follows from © 



that the reduced volume per particle, v 
be written as 



V/(Nr c ), 



v(t,p) = 1 — rftanh 



can 



(3) 



In order to understand the relation between water-like 
anomalies and phase transitions we investigate the be- 
havior of molecular volume while approaching the GSPT 
t = p = 0. Close to this point the value of the reduced 
volume depends on the path along which criticality is 
approached, namely 



lim lim v = 

p->0 t->0 

lim lim v = 1. 
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(4b) 



The behavior of the reduced volume indicates that 
slightly below and above the LDL-HDL critical pressure 
the system is arranged according the LDL and the HDL, 
respectively, while exactly at the critical pressure the sys- 
tem becomes a mixture of the two phases as the tempera- 
ture is reduced. Thus, the entropy gained by mixing with 
the HDL phase explains the density increase observed on 
the LDL side. 
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TABLE I. Coefficients of the core softened pair potential of 
Eq. ([5]). Distances in units of a = o\. 

Next, we consider if the presence of three critical points 
and two regions of density anomaly observed for d = 1 
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FIG. 4. Pressure versus temperature phase diagram for the 
3D case. The lighter lines are isochores for p* = 0.12 — 0.40, 
the solid lines the TMD lines and filled circles are the (esti- 
mated) location of the critical points. 



holds for d > 1. To this end, we perform molecular dy- 
namic simulations, using the HOOMD-blue packag e 26 ' 27 , 
at constant temperature (NVT) and constant pressure 
(NPT) in a three dimensional system composed of 2048 
particles in a cubic box with periodic boundary condi- 
tions, interacting through a continuous pair potential ob- 
tained by the addition of 3 different Fcrmi-Dirac distri- 
butions, 



U 



4 



A, 



exp 



(^) 



(5) 



with a — 0.08cr and the other coefficients given in TableQ] 
The potential illustrated as the smooth curve in Fig. Q] 
was chosen to mimic the three length scales and energies 
scales employed in the one dimensional case. The three 
length scales are identified by the change of the slope in 
the forced. The accessibility of the three length scales 
and the presence of the density anomalous region in the 
stable region of the pressure-temperature phase diagram 
is obtained by making the forces smooth and positive 2 ^. 
Units are reduced as in the ID case, except for pressure 
and particle density which are reduced as P* = Pa 3 /e 
and p* = a 3 p. 

The pressure versus temperature phase diagram is 
shown in Fig. |4j At low temperatures, the system 
presents gas, LDL, HDL, and VHDL phases. Gray lines 
arc the isochores p* = 0.12 — 0.40 and filled circles 
indicate gas-LDL, LDL-HDL, and HDL- VHDL critical 
points. Similarly to what happens in the ID case, in 
the vicinity of both LDL-HDL and HDL- VHDL critical 
points there are TMD lines located in the stable region of 
the pressure-temperature phase diagram. In both cases 
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the TMD lines seems to approach the critical pressures^, 
allowing us to test if the oscillatory behavior observed in 
a in the one dimensional case can also occur in the 3D 
case. This behavior is observed in Fig. [5l where the ther- 
mal expansion coefficient is drawn as a function of the 
pressure for different fixed temperatures T* — 0.48 — 0.68 
in the vicinity of the LDL-HDL critical point. As in the 
ID case a goes to — oo as P — >• P~ while it diverges to 
+00 as P — >• P c + at the critical temperature. 
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FIG. 5. Thermal expansion coefficient versus pressure for 
T* = 0.48 - 0.68. 



In this paper we used an exactly solvable one dimen- 
sional lattice model and molecular dynamic simulations 
of a three dimensional fluid model to show that three 
length scales core-softened potential can be designed to 
exhibit three liquid phases (LDL, HDL and VHDL), two 
liquid-liquid critical points and two density anomalous re- 
gions. In addition, we propose that the oscillatory behav- 
ior observed in the thermal expansion coefficient close to 
criticality can be considered as a signature of the connec- 
tion between a liquid-liquid critical point and a region of 
anomalous density increase, in the studied models. In the 
particular case of water the hypothesis of a liquid-liquid 
critical point is still under debate due to the impossibil- 
ity to experimentally probe the system in its expected 
location, at temperatures below the homogeneous nucle- 
ation temperature. We also propose that measurements 
of a as a function of pressure could be performed at tem- 
peratures well above criticality, where the system can be 
accessed (even though metastable), and used as a tool 
to test for the presence of a metastable second critical 
point. 
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